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A  code  was  developed  using  Newtons  iteration  applied  to  the  steady  Navier-Stokes 
equations.  Roe's  upwind  flux  difference  splitting  was  used  for  spatial 
discretization.  Second  order  accurate  values  of  the  Riemann  variables  were 
obtained  by  Spekreijses  interpolation  procedure  with  Van  Albada's  limiter  used 
to  reduce  spurious  oscillations  in  the  results.  The  Jacobian  derivatives  and 
the  code  necessary  to  compute  the  Jacobian  elements  were  generated  by  use  of 
the  symbolic  manipulation  code  MACSYMA.  The  Boeing  package  RSLIB  was  used  to 
invert  the  Jacobian  matrix. 

Computational  results  were  obtained  for  flat  plates,  flat  plate-ramp,  and  ogive- 
tangent-cvlinder  configurations  at  Mach  numbers  ranging  from  2  to  14.1.  These 
cases  contained  shock  on  shock  interactions,  viscous  layer  separations,  and 
turbulent  boundary  layer  transitions.  Convergence  was  obtained  in  15  to  200+ 
iterations.  The  results  agreed  well  with  experiment  or  other  computation  and 
demonstrated  that  Newton  iteration  can  produce  results  for  complex  high  speed 
flows  with  shock  waves.  • 
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Introduction 


This  document  comprises  the  final  technical  report  for  ARO  Grant  DAAL03-86-G-0039, 
"Graduate  Feiiowship  in  Hypersonics."  The  fellowship  was  awarded  to  Mr.  Paul  D.  Orkwis  after 
national  advertisement.  Mr.  Orkwis  completed  the  Requirements  for  the  Ph.D.  degree  on 
November  5,  1990  and  will  be  awarded  this  degree  at  NCSU's  December  19,  1990  graduation. 

He  has  accepted  a  position  as  Assistant  Professor  in  the  Department  of  Aerospace  Engineering  at 
the  University  of  Cincinnati  beginning  in  January  1991. 

The  remainder  of  this  document  will  discuss  the  research  problem  addressed  and  will 
present  a  brief  summary  of  research  results.  A  list  of  publications  and  presentations  of  the  work  is 
included. 

Research  Problem 

The  Ballistic  Research  Laboratory  (BRL)  at  the  Aberdeen  Proving  Ground  has  been 
studying  a  class  of  ballistic  "sustainer"  projectiles  which  seek  to  reduce  drag  by  igniting  gas- 
producing  solid  fuel  in  the  base  of  the  projectile.  In  cooperation  with  BRL,  developing  innovative 
means  for  computing  flow  over  these  projectiles  was  chosen  as  the  research  project  for  this 
fellowship.  BRL  personnel  have  kindly  supplied  advice  and  computational  resources  in  support  of 
the  project. 

These  projectiles  are  essentially  secant-ogive-cylinder- boattail  (SOCBT)  shapes^1 )  with  a 
bluff  base  region  where  the  combustion  products  are  injected  into  the  flow.  The  projectiles  travel 
at  supersonic  speeds  at  a  small  angle  of  attack  during  their  trajectory.  This  presents  a  difficult 
computational  problem  in  that  the  flow  is  three-dimensional  and  the  governing  equations  are  very 
stiff  due  to  the  source  terms  and  the  widely  varying  eigenvalues  in  the  free  stream  and  base 
regions.  Research  was  initially  directed  toward  finding  means  of  integrating  the  governing 
equations  accurately  and  efficiently  in  spite  of  this  stiffness. 

Research  Results 

Fishers  equation  provides  an  excellent  model  problem  for  testing  integration  techniques  for 
the  Navier-Stokes  equations  with  source  terms.  Fishers  equation  has  both  diffusive  and  source 


terms  and  has  an  exact  solution  for  comparison  with  the  computational  results.  A  combined 
integration  technique  was  applied  to  this  equation  to  allow  complete  freedom  to  vary  the  integration 
technique  continuously  between  fully  explicit  and  fully  implicit.  The  scheme  was  also  applied  such 
that  the  discretization  of  the  diffusion  term  and  source  term  could  be  different.  A  large  number  of 
discretization  combinations  were  tried  for  various  marching  step  and  source  term  magnitudes. 
Results  obtained  for  fully  explicit  discretization  were  accurate  for  small  source  terms  with  time 
steps  less  than  or  equal  to  the  stability  limit  but  were  of  course  unstable  for  larger  time  steps. 
Trapezoidal  (TR)  differencing  was  found  to  match  the  small  timestep  fully  explicit  results  but  with 
significantly  higher  timesteps.  Fully  implicit  (FI)  differencing  also  allowed  large  timesteps  but 
resulted  in  significant  phase  errors,  '"his  study  found  that  the  most  accurate  and  efficient  technique 
for  diffusive  equations  with  source  terms  was  the  TR  differencing. 

Studies  were  then  performed  with  the  TR  and  FI  schemes  on  more  Navier-Stokes-like 
model  equations.  Stability  and  modified  equation  analyses  of  the  linear  Burgers  equation  showed 
that  virtually  unlimited  timesteps  could  be  used  with  the  FI  scheme  but  not  with  the  TR 
differencing.  These  results  indicated  that  the  most  efficient  way  to  obtain  a  solution  of  given 
accuracy  of  the  linear  Burgers  equation  was  to  use  the  FI  scheme  with  as  large  a  timestep  as 
possible.  Convergence  in  as  little  as  2  iterations  was  obtained  with  the  FI  scheme. 

Tests  were  then  conducted  using  the  FI  scheme  for  the  nonlinear  Burgers  equation:  which 
indicated  that  reduced  iteration  counts  would  result  from  using  large  timesteps.  However, 
obtaining  a  given  level  of  convergence  could  not  be  guaranteed  for  the  nonlinear  equation  set.  It 
was  determined  that  this  difficulty  was  the  result  of  applying  a  solver  designed  for  linear  problems 
to  nonlinear  equation  sets.  Fortunately,  a  corresponding  solver  for  nonlinear  equation  sets  does 
exist:  Newton's  method.  This  method  was  then  used  for  the  next  steps  in  the  model  problem 
hierachy,  the  Euler  equations  and  finally  the  Navier-Stoke^  equations. 


The  Newtons  method  integrates  the  full  steady  Navier-Stokes  equations.  An  algebraic 
turbulence  model  is  included  to  simulate  fluid  turbulence  where  appropriate.  Newtons  method  as 
used  in  the  present  workU)  can  be  written 


3F(U) 

au 


AU  =  -  F(U) 


where  F(U)  is  the  set  of  discretized  Navier-Stokes  equations  and  the  vector  U  consists  of  the 


conserved  variable  vectors  at  each  point  in  the  discretized  domain.  The  Jacobian  Matrix 


3F(U)  . 
— I — 1S 

au 


taken  after  discretization  of  the  equations  and  is  exact  without  approximation  in  the  present  work. 
Roes  flux  difference  splitting  was  used  to  discretize  the  equations  since  an  analytical  expression  for 
the  spatial  discretization  was  necessary  to  obtain  an  exact  Jacobian.  Spekreijse's  interpolation 
procedure  was  used  to  obtain  the  left  and  right  states  for  Roe's  method  to  second  order  accuracy. 
Van  Albada's  limiter  was  used  to  prevent  non-physical  oscillations  in  the  solution. 

The  symbolic  manipulation  code  MACSYMA  was  used  to  perform  the  differentiations  and 
output  the  FORTRAN  code  necessary  to  compute  the  elements  of  the  Jacobian  matrix. 

MACSYMA  macros  were  written  that  produced  the  approximately  8000  lines  of  code  necessary  for 
each  differencing  stencil  in  approximately  one  hour  of  mini-computer  time.  The  use  of 
MACSYMA  allowed  very  rapid  formulation  changes  and  guaranteed  that  the  code  was  error  free. 

The  Newton  iteration  procedure  was  performed  by  successively  inverting  the  Jacobian 
matrix,  solving  for  AU  and  then  updating  Un+1  and  the  Jacobian.  No  attempts  were  made  to 
accelerate  convergence  due  to  complexity  of  the  compressible  flowfields  with  shock  waves.  The 
code  RSLIB  developed  at  Boeing  was  modified  to  be  used  with  the  ILU(0)  iterations  and  was 
used  to  perform  the  inversion  of  the  sparse  Jacobian  matrix. 

A  time  term  was  added  to  the  diagonal  of  Newtons  method  to  allow  initial  guesses  for  the 
solution  that  were  not  "close"  to  the  final  solution.  This  term  essentially  made  Newtons  behave 
like  a  backwards  Euler  formulation  that  would  approach  Newtons  method  as  the  time  term  was 
increased.  Slug  flow  initial  conditions  could  then  be  used  without  difficulty.  A  procedure  based 


on  the  norm  of  the  residual  was  developed  to  increase  At  automatically  during  the  iteration 
procedure.  Further  details  of  the  method  and  iteration  procedure  are  available  in  References  2.  3 
and  4. 

Computational  results  were  obtained  for  flat  plate,  flat  plate-ramp,  and  ogive-tangent- 
cylinder  (OTC)  configurations  at  Mach  numbers  from  2.0  to  14.1.  These  cases  contained  laminar 
and  turbulent  viscous  layers,  separated  regions,  shock  on  shock  interaction  and  strong  expansions. 
Convergence  to  a  residual  value  of  10'7  was  obtained  at  iteration  counts  ranging  from  16  to  217 
iterations.  A  discussion  of  two  of  these  cases  follows. 

A  flat  plate-ramp  configuration  tested  experimentally  by  Holden  and  Moselle  has  been 
used  often  for  code  verification.  Flow  at  ramp  angles  of  15°  and  18°  contain  separated  flow 
regions  (18°)  and  shock  on  shock  interactions  (15°  and  18°)  to  provide  a  severe  test  case  for 
Newton  iteration.  Results  were  obtained  for  the  15°  case  using  an  103  x  40  grid  comprised  of 
three  overlapping  40  x  40  grids  in  38,  62,  and  106  iterations  respectively.  Examination  of  the 
computational  results  revealed  excellent  agreement  with  the  experimental  pressure  distributions  and 
good  agreement  with  the  experimental  skin  friction  distribution  for  both  the  15°  and  18°  (which 
used  a  142  x  40  grid  and  required  38,  76,  162,  and  69  iterations  respectively)  cases. 

Flows  over  the  OTC  configuration  tested  experimentally  in  ref  (1)  were  also  computed. 
Five  overlapping  grids  of  40  x  40  were  used  for  an  overall  176  x  40  grid.  The  axisymmetric 
equations  proved  to  require  more  iterations  to  converge  the  solution  and  also  required  the  addition 
of  Harten’s  eigenvalue  correction  to  remove  difficulties  encountered  at  the  sharp  nose  of  the 
projectile.  'Hie  surface  pressure  data  agreed  well  between  experiment  and  computation  except  near 
transition.  Further  work  is  needed  to  investigate  relaxation  models  in  this  region. 

Concluding  Remarks 

This  work  provides  the  initial  demonstration  that  Newton  iteration  can  be  used  to  obtain 
computational  results  for  complex  compressible  flows  including  strong  shock  waves.  It  has  also 
been  demonstrated  that  the  code  to  compute  exact  Jacobians  can  be  easily  obtained  through  use  of 
MACSYMA  to  allow  rapid  generation  of  code  for  new  formulations.  Newton  iteration  is  not  yet 


competitive  in  efficiency  with  standard  time  marching  techniques.  However,  future  research 
should  improve  efficiency  by  use  of  convergence  acceleration  and  improved  sparse  matrix 
inversion  packages. 
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